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1 .  Introduction 


From  classical  potential  theory  we  know  that  the 
(Newtonian)  potential  of  the  earth's  body  E 

V  (x)  =  4  J  p  (y)  dv(y)  (l.D 

E 

(  4.  :  gravitational  constant) 

is  a  harmonic  function  in  free  space.  That  is,  V  is 
arbitrarily  often  differentiable  in  the  earth's  exterior 
Ee  and  it  satisfies  Laplace's  equation 

AV(x)  *  0  (1.2) 

for  all  points  x  €  e^.  Moreover,  V  is  regular  at 
infinity,  i.e.  V  and  VV  tend  to  zero 

|V(x)  I  *  0(j^-) 

1 V  V  (x)  I  «  o  <-*-*> 

Ixl2 

if  | x |  tends  to  infinity.  Since  the  density  function  p  is 
better  known  only  for  parts  of  the  upper  crust  of  the 
earth,  the  closed  representation  (1.1)  of  the  potential  V, 
however,  cannot  be  used  in  numerical  computations.  Instead, 
we  have  to  look  for  suitable  approximations.  In  this 
connection  a  result  first  formulated  by  C.  Runge  (1885) 
and  later  generalized  by  J.L.  Walsh  (1929)  and  T.  Krarup 
(1969)  is  of  particular  importance.  It  states: 


(1.3) 
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2 


Any  function  V  harmonic  outside  the  earth's 
surface  S  and. regular  at  infinity,  may  be 
approximated  by  a  function  U,  harmonic  outside  an 
arbitrarily  given  sphere  (Bjerhammar  sphere)  inside 
the  earth  and  regular  at  infinity  in  the  sense  that 
for  any  given  e  >  0,  the  relations 

|V{x)  -  U(x)  j  <  E 

resp.  (1’4) 

|  VV(x)  -  VU(x)  |  <  e 

hold  for  all  points  x  of  the  Euclidean  space 
outside  and  on  any  closed  surface  completely 
surrounding  the  earth's  surface  in  the  outer  space. 

The  value  e  may  be  arbitrarily  small,  and  the 
surrounding  surface  may  be  arbitrarily  close  to  the 
earth's  surface  (cf.  Moritz  (1980)). 

In  this  formulation  the  Runge-Walsh-Krarup  theorem  is 
a  pure  existence  theorem,  i.e.  it  guarentees  only  the 
existence  of  an  approximating  function  U  and  does  not 
provide  a  method  to  find  it.  Nothing  is  said  about  the 
structure  of  the  approximation  U.  The  theorem  describes 
merely  the  theoretical  background  for  practical  purposes 
of  determining  an  approximation  to  the  external 
gravitational  potential  of  the  earth. 

In  a  spherically  symmetric  model,  however,  the 
situation  is  completely  different.  For  a  spherical  earth, 
constructive  approximations  of  V  are  available  by  use  of 
trial  functions  with  a  larger  harmonicity  domain.  As 
trial  functions  the  following  (base)  systems  {<J>  }  can  be 
used,  for  example: 


solid  spherical  harmonics 


IS  }  maximal  linearly  independent  system  of 

surface  spherical  harmonics  (origin  at  earth's 
center  of  gravity) 

(ii)  "  mass  points  " 


{Xp}  "fundamental"  sequence  of  points  in  the 
complement  of  Eg  U  S  (of.  Chapter  8) 


(iii)  "  multipoles  " 


,  3  \  1 

3V 


[  j]  *n 


xQ  is  a  fixed  point  in  the  complement  of  Eg  u  S 
j  =  (^^2^3)  is  a  triple  of  non-negative 
integers , 


t  jl  j 1 2+j 3  •'  (3x  )  = 


ax13l3x2:j29x3:,3 


We  recapitulate  the  standard  method  of  approximating  an 
external  potential  V  in  a  spherical  model. 


Let  us  suppose  that  S  is  a  sphere  with  radius  R  about 
the  origin  and  E  is  the  outer  space  of  S: 

C 

S  =  (x  6  R3  i  i x |  —  R} 

Ee  =  {x  €  R3  [  |x|  >  R}. 

Let  V  be  continuous  in  Eg  U  S  and  twice  continuously 

differentiable  in  E  .  Moreover,  let  V  satisfy  AV(x)  =  0 

e 

for  all  x  €  E  . 

e 

Let  us  denote  by  {0^}  one  of  the  linearly  independent 
systems  of  trial  functions  listed  above.  Corresponding  to 
the  countably  infinite  sequence  {0n>  there  exists  a  system 
{$*}  of  trial  functions  <f>*  orthonormalized  with  respect  to 
the  inner  product  of  the  space  *2(S)  of  all  square 
integrable  functions  on  S  (Gram-Schmidt  orthonormalizing 
process) 


/  <t>*  (y)ds  (y)  =  5^  . 

s 


(dS  i  surface  element) 


(1.5) 


Then  the  potential  V  can  be  represented  by  the  expansion 


V(x)  =  l  (V , <t>*)  o*(x)  , 
n=o  n 


(1.6) 


where  the  numbers 


(V,<j>*)  =  /  V  (y)  4>*(y)dS(y) 
S 


(1.7) 


are  the  Fourier  (or  orthogonal)  coefficients  of  V  on  S 
with  respect  to  the  '  -stem  f  ■*  } . 
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More  explicitly  this  reads:  given  an  error  bound  e  >  0, 
then  there  exists  an  integer  N  =  N(e)  such  that 

N 

|v{x)  -  l  (V,<J>*)  0*(x)  |  <  s  (1.8) 

n=o 

for  all  points  x  €  G  with  G  <=  E  and  positive  distance  to  S 
(dist(G,S)  >  A  >  0) .  In  each  compact  subset  G  c  Ee  the 
convergence  is  uniform. 

Therefore  we  have  a  constructive  version  of  Runge ' s 
theorem  by  means  of  the  truncations 

,Nv  N 

VKNI  (x)  =  l  (V,4>*)  <D*(x)  (1.9) 

n=o 

in  a  spherical  (earth's)  model. 


The  approximation  technique  is  formulated  so  as  to 
have  the  permanence  property:  the  transition  from  v(N) 
t /~  V(N+1)  necessitates  merely  the  addition  of  one  more 
term,  all  the  terms  obtained  formerly  remaining  unchanged. 


Furthermore 

/  I V  (y) 

S 


(N) 

V  has  the  following  minimum  property 


N 

-  I  (V,d>*)  4>*(y) 
n=o 


dS 


(l.lo) 


I 

S 


V  (y) 


“  I  an  0* (y) 1 2  ds 

n=o 


for  any  selection  of  coefficients  an> 


The  considerations  given  above  lead  to  the  impression 
that  a  least  square  approximation  by  truncated  expansions 
V(N)  of  the  form  (1.9)  is  intimately  related  to  spherical 
boundaries.  Our  purpose,  however,  is  to  prove  that  the 
essential  steps  leading  to  a  least  square  approximation 
of  an  external  gravitational  potential  by  truncations  of 
series  expansions  into  (multi) poles  can  be  transcribed 
in  the  same  way  to  general  (non-spherical)  models. 

We  summarize  the  main  results. 


Not-necessarily  spherical  models : 

Let  us  suppose  that  S  is  the  boundary-surface  of 
a  not-necessarily  spherical  earth's  model  (e.g. 
ellipsoid,  telluroid,  geo  id ,.  real  (regular)  earth's  surface) . 
Let  Ee  be  the  outer  space. 

Let  V  be  continuous  in  E  US  and  twice  continuously 

e 

differentiable  in  Eq.  Moreover,  let  V  satisfy 
AV(x)  =  0  for  all  x  e  Eg. 

Let  us  denote  by  {<fcn}  one  of  the  linearly  indepen¬ 
dent  systems  listed  above.  Corresponding  to  the 
sequence  { n }  there  exists  a  system  of 

trial  functions  <p*  orthonormalized  with  respect 

n  2 

to  the  inner  product  of  the  space  (S)  of  all  square 
integrable  functions  on  S.  Then  the  potential  V  can 
be  respresented  by  the  expansion 

oo 

V (x)  =  l  (V, 4> *)  0*(x)  , 

n=o 


(1.11) 


where  the  numbers 


(V,<D*)  =  /  V  (y)  <D*(y)  dS(y)  (1.12) 

S 

are  the  Fourier  (or  orthogonal)  coefficients  of  V  on  S 
with  respect  to  the  system  {<$>*}. 

More  explicitly:  given  an  error  bound  e  >  0,  then 
there  exists  an  integer  N=N ( e )  such  that 

N 

I V (x)  -  l  (V,**)  <J»*(x)  |  <  e  (1.13) 

n=o 

for  all  x  €  G  with  G  c  Eg  and  dist(G,S)  >  X  >  o.  In 
each  compact  subset  G  c  Ee,the  convergence  is  uniform. 


It  should  be  pointed  out 
that  the  series  (1.11) 
guarentees  ordinary  point- 
wise  convergence  not  only 
outside  the  smallest  bounding 
sphere,  but  indeed  in  the 
(whole)  outer  space  Ee- 
Our  approach  therefore  avoids  inadequate  phenomena  of 
divergence  in  regions  between  a  bounding  sphere  and  the 
surface  S.  The  price  to  be  paid  is  the  orthonormalization 
process  with  respect  to  the  not-necessarily  spherical 
surface  S. 

As  in  the  spherical  case,  the  advantages  are  the 
permanence  property  and  the  minimum  property. 


At  the  surface  S,  the  approximation  by  trial  functions 
depends  closely  on  the  smoothness  presupposed  for  V. 

Provided  that  V  is  only  assumed  to  be  continuous  on 
the  surface  S,  then  it  can  be  proved  that,  to  every  e  >  0, 
there  exist  an  integer  N  =  N(e)  and  coefficients  b  , ...,bN 
such  that 

N 

sud  'V(x)  -  l  b  4>*(x)|  <  e  (1.14) 

x€S  n=o 

and  by  maximum  principle  of  potential  theory 

N 

sup  |V(x)  -  l  b  <j>*(x)  |  v  e  (1.15) 

xGE  uS  n=o 

e 

(cf.  Freeden ( 1 980)  Chapter  6).  But  the  coefficients  b0,...,fc^ 
do  not  coincide,  in  general,  with  the  Fourier  coefficients 
(V , <p* ),...,  (V,4>jUp  . 

Hence,  "downward  continuation"  of  the  expansion 

oo 

I  (V,**)**  (1.16) 

n=o 

to  the  surface  S  cannot  be  assured  generally  under  the 
assumption  of  continuity  of  V  over  the  surface  S. 

The  series  expansion 

o° 

V (x)  =  l  (V, $*)$* (x) ,  (1.17) 

n=o 

however,  remains  true  for  all  x  €  Ee  U  S  in  the  ordinary 
pointwise  sense  if  stronger  conditions  of  smoothness  are 
imposed  on  V  in  addition  .  In  the  spherical  case,  for 
example,  the  continuous  dif f erentiahilitv  of  v  on  .c  is 
sufficient  to  verify  the  validity  cf  (1.17)  on  c  for 
spherical  harmonics.  (  detailed  discussion  of  "^o*’n.'-’ard 
continuation"  of  series  expansions  into  trial  furcti cns 
for  non-snherical  cases  seems  to  be  oner) . 


In  detail  we  are  concerned  with  the  following 
investigations : 


In  Chapter  2  and  3  notations  are  introduced. 

Essential  tools  of  the  presentation  are  the  "jump 
relations"  and  the  "limit  formulas"  of  potential  theory 
(Chapt.4)  and  a  regularization  theorem  (Chapt.5).  The 
adventages  of  an  approximation  by  harmonic  trial  functions 
are  discussed  in  Chapt.6  and  7.  Then  non-constructive 
approximation  theorems  (Chapt.8)  and  constructive 
approximation  theorems  (Chapt.9)  are  developed.  Chapt.10 
gives  an  application  to  the  external  gravitational 
potential  of  the  earth.  In  Chapt.11  a  numerical  method 
based  on  the  Cholesky  factorization  is  presented.  The 
efficiency  of  the  numerical  method  is  analyzed.  Finally 
( Chap t. 1 2) ,  some  simple  examples  give  an  impression  of 
the  accuracy. 
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2.  Basic  Settings  of  Functional  Analysis 

The  following  definitions  and  notations  are  standard 
in  functional  analysis. 

A  linear  space  is  a  set  of  elements  such  that  to  any 

finite  subset  the  linear  combination 

o  n 

a  (p  +  . . .  +  a  <p  (a  :  scalars)  is  also  a  member  of  the 
o  o  n  n  n 

linear  space.  We  shall  only  consider  real  linear  spaces. 
Hence  aQ/...  ,  an  are  assumed  to  be  real.  Any  linear 
space  possesses  a  unique  zero  element .  A  collection 

of  elements  out  of  a  linear  space  is  called 
linearly  dependent  if  there  exist  numbers  a0,...,a  ,  not 
all  of  them  being  zero,  such  that 

a  +  .  .  .  +  a  <p_  =  0 .  (2.1) 

o  o  n  n 

Otherwise  the  elements  <PQ,...,(Pn  are  called  linearly 
independent . 

For  each  subset  of  a  linear  space  we  consider  the 
space  generated  by  this  subset,  i.e.  linear  combinations 
of  the  form 


Vo  + 


+  °n*n 


The  union  of  all  such  spaces  is  called  the  span  of 


A  normed  space  X  is  a  linear  space  which  assigns  to 
any  of  its  elements  (p  a  non-negative  number  II  (pll,  called 
the  norm  of  (p,  with  the  following  properties: 


(i)  II  (pll  >  0;  II  (pll  ■  0  «  (p  =  0 

(ii)  II  tp  +  (fill  <  II  (pll  +  II  (fill 

(iii)  II  a(pll  =  |  a  |  II  <pll  . 


(a :  real)  . 


The  norm  describes  distances  between  pairs  of  elements 
(P , vp  £  X.  The  distance  between  <P  and  <P  is  defined  by 


dist  (<P  ,i fi  )  =  !l(p  -(fi  I 


(2.2) 


The  distance  between  two  subsets  A,B  of  X  is  defined  by 


dist(A,B)  =  inf  II tp  -  tfill  . 
(p€A 
iptB 


(2.3) 


Distances  can  be  used  to  define  convergence.  Two  cases 
are  of  particular  interest. 

(1)  Convergence  to  a  limit  element . 

Let  {(P  }  be  a  sequence  of  elements  out  of 

“  n=o, 1 , . . . 

a  normed  linear  space  X.  Let  (p  be  another  element 
of  X.  We  say  that  the  sequence 

n  n=o , 1 _ _ 

converges  to  the  limit  element  (p,  in  symbols. 


lim  =  <p 
n-«  n 


or  (Pn  -*■  ^  in  X  , 


if,  given  e  >  0,  there  exists  an  integer  N  =  N(e) 
such  that 
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to  -  tpll  <  e 
n 


(2.4) 


for  all  n  >  N(e) 


(ii)  Convergence  in  the  sense  of  Cauchy 

A  sequence  { tp  }  is  called  convergent  in  the 

n  n=o , 1 , . . . 

sense  of  Cauchy,  if  given  e  >0,  there  exists  an 
integer  N  *  N(e)  such  that 


%  -  v  4  c 


(2.5) 


for  all  n,m  >  N (e) . 

We  see  that  a  Cauchy  sequence  is  defined  without  using  a 
limit  element.  It  merely  requires  that  the  elements  of 
the  sequence  cluster  in  a  certain  way. 

A  normed  space  X  is  called  complete,  if,  to  every 

Cauchy  convergent  sequence  {ip„}  in  X,  there 

n  n=o , 1 , . . . 

exists  an  element  ip  €  X,  such  that  {ip  }  converges 

n  n=o , 1 , . . . 

to  the  limit  element  <p. 

As  soon  as  the  elements  cluster  in  the  sense  of  Cauchy 
in  a  complete  normed  space,  a  limit  element  must  exist. 

Complete  normed  spaces  are  called  Banach  spaces . 


For  a  normed  space  X,  a  system  (cp  }  is  called 

n  n=o ,  1  ,  .  .  . 

a  basis ,  if  for  any  e  >  0,  ip  €.  X  there  exist  an  integer 
N=N(e)  and  coefficients  a  a.,  such  that 

O  ,  .  .  . ,  N 


tp  -  l  oij^ll  <  e. 
k=o 


if  to  any  pair  of  its  members  a  real  number  , 

called  the  inner  product,  is  assigned  which  has  the 
following  properties: 


(i) 

(<P,  ip) 

(ii) 

(cup, tp)  = 

a 

(iii) 

(<P1  +q>2 ,  ^ ) 

=  (ip1f\p)  +  (ip2 

(iv) 

(ip,<p)  >  0 

,  (ip,tp)  =  0  ♦* 

real) 


A  consequence  of  these  properties  is  the  Schwarz 
inequality 

(<P,<i>)2  <  (<p,<p)  (ip  ,<p)  .  (2.6 

We  say  that  a  system  {<p.  }  „  of  functions  of  an 

^  k — o ,  1  ,  • »  . 

inner  product  space  X  is  an  orthogonal  system  if  its 
elements  are  mutually  perpendicular,  that  is 

(^k'^i*  =  0  (2.7 

if  i  *  k,  and  in  addition 

(^k'^k*  *  0 


for  all  k. 

An  orthogonal  system  {ip  }  is  called 

n  n=o, 1 , . . . 

orthonormal  system  if  (<pk,<p^)  =  1  for  all  k.  One  can 
always  obtain  an  orthonormal  system  by  dividing  each 
element  by  its  norm. 


A  further  property  of  the  inner  product  which  may  be 
deduced  from  the  defining  properties  (i)  -  (iv)  is  the 
continuity  of  the  inner  product .  This  means  that  if 
<Pn  ♦  <Pm  ■*  ipf  then  (<Pn,(J>m)  -*■  (<p, ip). 

The  inner  product  induces  a  norm  by  defining 

11  (pH  »  ✓  ((p'#cp~)  .  (2.8) 

Hence,  an  inner  product  space  X  is  automatically  a 
normed  space.  Schwarz's  inequality  may  be  reformulated 
by  taking  the  square  root 

|  I  <  II  tpll  *  lityll  .  (2.9) 

A  Hilbert  space  is  an  inner  product  space  which,  in 
addition,  is  complete. 

If  a  basis  {<p„)  in  an  Hilbert  space  X  is 

n  n=o, 1 , . . . 

orthonormalized,  each  element  ip  of  the  space  X  is 
representable  by  the  generalized  Fourier  series 

oo 

‘ o-l  (<P^n)<Pn  ( 2 . 1o) 

n=o 

in  the  sense  of  X. 

An  orthonormal  basis  in  a  Hilbert  space  is  often 
called  Hilbert  basis. 


3.  The  Banach  Space  (S)  and  the  Hilbert  Space  (S) 

3 

If  G  is  a  set  of  points  in  R  ,  3G  will  represent  its 
boundary.  The  set  G  =  G  U  3G  will  be  called  the  closure 
of  G. 

A  set  G  c  R3  is  called  a  region  if  it  is  open  and 
connected . 

The  restriction  of  a  function  f  to  a  subset  M  of  its 
domain  is  denoted  by  f  jM,  for  a  set  I  of  functions  we 
define  l|M  as  the  set  of  restrictions  f|M: 

I | M  *  {f | M| f  €  1} 

A  function  f  possessing  k  continuous  derivatives  is 
said  to  be  of  class  C  . 

A  surface  S  c  R3  will  be  called  regular,  if  it 
satisfies  the  following  properties: 

3 

(i)  S  divides  three-dimensional  Euclidean  space  R 
into  the  bounded  region  E  *=  Ei  (inner  space)  and 
the  unbounded  region  E  (outer  space)  defined  by 

3  _  e 

E@  -  R  -  E 

(ii)  S  is  a  closed  and  compact  surface  with  no  double 
points 

(iii)  S  is  a  -surface 

From  the  definition  it  is  clear  that  all  (geodetically 
relevant)  earth's  models  are  included.  Regular  surfaces 
are  for  example  sphere,  ellipsoid,  spheroid,  telluroid, 
geoid  or  real  (regular)  earth's  surface. 
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4.  Limit  Relations  and  Jump  Relations  in  Potential 
Theory 

At  each  point  y  of  the  boundary  S  =  3E  we  can  construct 
a  normal  n(y)  pointing  into  the  outer  space  Eg.  The  set 

ST  =  {x  6  R3  |  x  =  y  +  tn (y) ,  y 6  S>  (4.1) 

generates  a  parallel  surface  which  is  exterior  to  E  for 
t  >  0  and  interior  for  t  <  0.  It  is  well-known  from 
differential  geometry  (cf.  Cl.Mtiller  ( 1  969)  ,§  13)  that 
if  | r |  is  sufficiently  small,  then  the  surface  S  is 
regular,  and  the  normal  to  one  parallel  surface  is  a 
normal  to  the  other. 

Let  g  be  a  continuous  function  on  S  (g  e  (S) ) . 

Then  the  functions  UN  defined  by 

UH<*>  -  1  q(y)(^j|-r)N  — i-r  dS(y)  (4.2) 

o 

are  infinitely  often  differentiable  in  Ei  and  Eg  and 
satisfy  the  Laplace  equation.  In  addition  r  they  are 
regular  at  infinity. 

The  function  defined  by 

^  (x)  =  J  g(y)  ]x-y\  ds^y) 

s 

is  called  the  potential  of  the  single  layer  on  S,  while 
U2  given  by 

°2(*>  -  /  s(y)  mv  w  dS(y) 

is  called  the  potential  of  the  double  layer  on  S. 
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For  g  €  C  ^  (S)  the  functions  UN  can  be  continued 
continuously  to  the  surface  S,  but  the  limits  depend 
from  which  parallel  surface  (inner  or  outer)  the  points  x 
tend  to  S . 

On  the  other  hand,  the  functions  UN  also  are  defined 
on  the  surface  S,  i.e.  the  integrals 

UN(X)  *  /  j^T-dSly) 

exist  for  x  €  S. 

Furthermore,  the  integral 

o;<x>  -  j  g(y>  ds|y)  (4-3) 

exists  for  x  €  S. 

However,  the  integrals  do  not  coincide  in  general  with 
the  inner  or  outer  limits  of  the  potentials. 

From  classical  potential  theory  (cf.  e.g.  O.D.Kelogg 
(1929),  J-  Schauder  (1931))  it  is  known  that  for  x€  S 

lim  U1  (x  ±  tn(x))  =  U1 (x)  (4.4a) 

t-*-o 

T  >0 

3U1 

lim  - — (x  ±  tn  (x) ) 

_  .  _  3n 
t-*-o 

t  >0 

lim  U2(x  ±  xn(x))  =  ±  2Trg(x)  +  U2(x)  (4.4c) 

t-*-o 
t  >0 


=  + 


2irg  (x)  +  Uj  (x) 


(4.4b) 


("  limit  relations  ") 


2o 


i 

[■ 


k 


M 

3 


i 


b 


and 


lim(U.j  (x  +xn(x))  -  (x  -  xn(x))  =  0  (4.5a) 

x-*-o 
x  >o 


3U. 

lim(-r— -!-(x  +  xn(x) ) 
x+o 
x  >o 


3  0. 
3rT 


(x 


xn(xS)  =-4irg(x)  (4.5b) 


lim(U2(x  +  xn(x)  )  -  U2  (x  -  xn(x)  )  )  =  4rrg(x)  (4.5c) 

x--o 

x>o 


3U2 

lim(3— ^(x  +  xn(x) ) 

x-o  9n 

x>o 


xn(x) ))  =  O. 


( 4 . 5d) 


("  jump  relations  ”) 


In  addition  ,  O.D.  Kellogg  has  proved  that  the  above 
relations  hold  uniformly  with  respect  to  all  x  e  S. 

This  means  that 

lim  sup  jU^fx  ±  xn(x))  -  U1 (x) |  =  0  (4.6a) 

x-*-o  xes 

x>o 

3U. 

lim  sup  |-»—!-(x  ±  xn(x) )  ±  2-ng  (x)  -Ui(x)|=0  (4.6b) 
x-*o  x€S  9 
x  >o 


and 


v 

k 

k 

r 
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lim  sup  |U1(x  +  xn(x)  -  U, (x  -  xn(x))|  =  O  (4.7a) 

T"*o  xes 

T>0 


au1  3u. 

lim  sup  I  J— -(x  +  xn (x) )  -  3—Mx  -  in  (x))  +  4Trg  (x)  i=C  4.7b) 
x-*-o  X€S  d  8n 

T  >0 


lim  sup  |  U  2  ( x  +  in  (x)  )  -  U,  (x  -  xn(x))  -  4irg  (x)|=0  (4.7c) 

t-o  xes  * 

t>o 


au2  3u2 

iim  sup  j  r--(x  +  in  (x)  )  -  -^(x  -  xn(x)  )|  =0.  (4.7d) 

T+o  xes  d  dn  ' 

t  >o 


Herein 


± 


xn(x)> 


n(x)  *  (7U)  (x  ±  m(x) )  . 


These  relations  can  be  generalized  to  the  Hilbert 
2 

space  ^  (S)  using  exclusively  functional  analytic 
means  (cf.  Freeden ( 1 98o) ,  Kersten  ( 1 98o) ) . 

For  square-integrable  functions  g,  i.e.  g  e/2(S),  we 
then  have 

lim  J  1 U1  (x  ±  xn(x)l  -  U^x)]2  dS(x)  =  0  (4.8a) 

x-*-o  S 

T>0 


lim  /  |-r—  (x  ±  x  n  ( x  )J  -  U,'  (x)  ±  2  tt  g  (x)|2  dS(x)=0  (4.8b) 

x-*o  S  dn 

x>o 


lim  /  jU2(x  ±  rn(x)j  -  U2  (x)  +2irg  (x)  j  2dS  (x)  =  0  (4.8c) 

T-K>  S 

T>0 


and 

lim  /  |U,  (x  +  xn(x))  -  U.  (x  -  tn  (x))j  2dS  (x)  =  0  (4.9a) 

t-o  S  1 

T>0 

9U1  9U- 

lint  /  laTT(x  +  Tn<x))  "  tt-(x  ~  fn(x))+47rg(x)  |  dS(x)=0  (4.9b) 
X+O  S  dn  dTl 

X  vO 


lim  /  |U2(x  +  xn(x))  -  U,(x  -  xn  (x)  )  -  4ng(x)  1 2dS(x)=0  (4.9c) 

x-*o  S 

x>o 


9U2  9U2 

lim  /  |g^-(x  +  xn (x)J  -  (x  -  xn(x))  |  <iS(x)  =  o.  (4.9d) 
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5.  Regularization  Theorem 

Let  ^(Eg)  be  the  space  of  functions  V  with  the 
following  properties: 

(i)  V  is  twice  continuously  differentiable  in  Ee  and 
continuous  in  Ee=  E«U  S,  i.e.  V£  (Ee)H  (Ee) 

(ii)  V  satisfies  Laplace's  equation  in  Ee  : 

AV (x)  =  0  for  x  £  Ee  • 

(iii)  V  is  regular  at  infinity: 

|v(x)|  =  ot-^i 

j  7V  (x)  |  =  0(— !— 5) 

ixl2 

if  I  X  I  -►  00  . 


Suppose  now  that  {0  }  is  a  sequence  of 

n  n=o , 1 , . . . 

functions  <j>n  belonging  to  the  class  ^(Ee)  »  We  denote  by 

{tp  }  the  system  of  restrictions  <p  =0  |S,n=o,1 , . . .  . 

n  n«o , 1 , . . . 


Assume  that  the  system  {tp  >  converges 

n  n=o , 1 , . . . 

uniformly  to  the  function  <p  on  S,  i.e.  for  given  e  >  O, 
there  exists  an  integer  N  =  N(e)  such  that 


sup  j<p  (x)  -  tp(x)  |  <  e  (5.1) 

X€S  n 

for  all  n  >  N.  Then  the  maximum  principle  of  potential 

theory  assures  that  the  sequence  { 0  >  converges 

n  n=o , 1 , . .  . 

uniformly  to  the  (unique)  solution  0  of  the  exterior 
(Dirichlet)  problem  0  6  £(Ee)  ,  0|S  =  <p. 
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In  other  words,  this  means  that 


slip  |  <P  (x)  -  <f>  (x)  j  <  e 
x€E 


(5.2) 


for  all  n  >  N. 

If  the  system  { ip  >  is  only  convergent  to  (p 

n  n=o , 1 , . . . 

in  quadratic  sense,  i.e.  for  given  e  >  0,  there  exists  an 
integer  N  =  N(e)  such  that 


/  |<p  (x)  -  ip  (x )  |  aS  <  e  , 
S  n 


(5.3) 


then  we  cannot  expect  in  general  that  the  sequence 
is  uniformly  convergent  in  the  whole 

n=o , 1 , . . . 

space  Ee=  EeU  S.  Our  intend,  however,  is  to  the  show  that 

the  quadratic  convergence  of  a  sequence  {y  }  to 

n  n=o , 1 , . . . 

a  function  tp  on  the  surface  S  implies  ordinary  pointwise 


convergence  of  the  sequence  {$  } 


n  . 

n  o ,  1  , . .  . 


to  the  solution 


<J)  of  the  exterior  (Dirichlet)  problem  for  each  point  x  of 

the  outer  space  E&.  In  addition  ,  for  each  compact 

subset  B  c  Ee  with  positive  distance  to  the  surface  S 

(dist  (B , S)  >  A  >  o)  the  convergence  of  (<f>  } 

.  n  n=o , 1 , . . . 

is  uniform. 


For  that  purpose  we  need  the  following  regularization 
theorem : 


Let  B  be  a  compact  subset  of  the  outer  space 
Ee  with  dist  (B ,S)  >  A  >  0. 

Suppose  that  4> , f  are  functions  of  the  class  ^(E  ) 

e 

with  4>|S  =  <p  ,  V  j  S  =  il>.  Then  there  exists  a  positive 
constant  C  =  C(S,B)  such  that 


sup  I  $  (X)  -  Y  (X)  '  <  C  Jj  ;  ip  (y  >  -  <My )  !  dS, 

x€B  v  s 


(5.4) 
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Proof .  By  the  smoothness  condition  imposed  on  S  =  9E, 
there  exists  the  classical  Green  function  for  the 
scalar  Dirichlet  problem  related  to  Eg)  i.e.: 


4>{x)  *  /  ^(y)-g^Jy)  ^  (x,y)ds  (y) 

s 

'F(x)  *  /  ii-(y)— -L-}  ^  (x,y)dS(y) 


for  all  x  €  E, 


Consequently,  we  have 


$(x)  -  V  (x)  »  /[<p(y)  -  1  37^)  ^  (x,y)dS(y) 

S 


for  all  x  €  E  - 
e 


C  =  C  (S  ,B)  * 


sup  [  /  f  (x,y)  |2dS(y)l 

x€B  S 


Then,  for  each  x  €  B,  we  have  by  the  Cauchy-Schwarz 
inequality 


<Mx)  -  'Hx)  |2  =  i  }I  4?(y?-^(y)l  f^~yj  £  (x,y)dS(y)| 

S 

<  J|<p(y)-<My)  |2dS(y)/|I|T-)^(x,y)|2dS(y) 


<  C2  /  !<p{ y)  -  ip(y)  !  2ds  . 
s 
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Thus 


SU£  |«(x)  -  'F(x)!2  <  C2  /  !  <p  (y)  -  <P  (y )  |  2dS  (y) .  □ 

xeB  s 


V  ' 


As  consequence  of  the  regularity  theorem  we  obtain 
the  following  Corollaries: 


Corollary  1:  Let  {<£  }  denote  a  sequence  in 

- —  n  n=o , 1 , . . . 

fp (E&)  •  Denote  by  { tp  }  the  restrictions 

n  n=o , 1 , . . . 

-  4>  I S  of  the  functions  <j>  to  the  surface  S. 
n  n  n 

Suppose  that  <j>  is  the  (unique)  solution  of  the 
Dirichlet  problem 

<J>  €  fi(Ee)  ,  <j>jS  »  cp  . 

If  now  ipn  -*■  <p  in  £  z  (S)  : 

lim  |tp  (x)  -  <p(x)  1 2dS  (x)  *  0, 
n-*»  V  S 

then  <t>n  +  <J>  uniformly  on  each  compact  subset  B  <=  E0 
with  dist (B,S)  >  X  >  0, 

lim  sup  !<t>  (x)  —  (x)  i  =0. 

n-*°°  x€B  n 


Proof .  According  to  the  regularization  theorem  we  have 


sup  |  <J>  (x)  -  <f>  _(x)  |  <  C  //  l\p  (y)  -  «o  (y)  j2dS(y) 
x€B  n  /  S  n 


for  n  =  o, 1 , . . .  . 


27 


Thus  the  £  -convergence 

^/fl<p(y)  -  «Pn(y)l2  ds(y)  -  o 

implies 

SU£  |<j>  (X)  -  <p  (X)  (  -  0  , 

x€B  n 

as  n  tends  to  infinity.  O 

Qc^-Ollary  2;  Let  {'4)*)  be  an  orthonormal 

2  nn=o,if... 

system  in  *T*(S),ip*  e  C[0)  (S) .  Denote  by  {<»*}  ,  the 

sequence  of  functions  <)>*€  #Ee)with  <t>*  |  S  =  ip  *  for  n=o , 1 , .  . . 

Suppose  that  4>  is  the  (unique)  solution  of  the 
(Dirichlet)  problem 

4>  €  ^(Ee)  ,  4)  |  S  =  q>  . 

Then  the  sum 

(N)  N 

#  (x)  *  l 

n~o 

N 

j- 

=  l  /  <P(y)«p*(y)dS(y)  <Js*(x) 
n=o  S  n  J  n 

converges  uniformly  to  <Hx)  as  N  -»  •  on  each  compact 
subset  B  c  Ee  with  dist(B,S)  >  X  >  o, 


lim  su£  I  <p  (x)  -  <J>  (x)  | 

N-*-°°  x€B 


=  lim  sup 
x€B 


N  r 

l  /  <p (y) <p* (y)ds (y) 
n=o  S  " 


<p  *  (x) 
'  n 


<P  (x) 


=  0. 


Proof .  According  to  the  regularization  theorem  it 
follows  that 


*  C 


A 


N 

z 

n=o 


J 

L  S 


<P(y)4>*  (y)dS  (y) 


<p*  (x) 


-  ip(x)  |  dS  (x) 


2 

As  element  of  the  Hilbert  space  *  (S)  ,  the  function  ip  can 
be  represented  by  its  Fourier  series  with  respect  to  the 
system  (cp* }  : 


lim 

N-*-* 


i  <P(y)<P*(y)ds(y) 


<p*  (x) 


~ <p(x )  |  2dS  (x) 


This  proves  Corollary  2. 


O 
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Before  giving  a  list  of  the  most  important  systems 
of  trial  functions  satisfying  the  assumptions  of  the 
Corollaries  stated  above,  we  recapitulate  some 
functional  analytic  aspects  concerned  with  the  closure 
and  completeness  of  sequences  in  (general)  inner  product 
spaces . 


3o 


A  sequence  {tpn>  is  called  closed  in  a  normed 

n=o ; 1 i » ■ • 

space  X  with  norm  II  *11  ,  if,  given  ip  €  X  and  e  >  0,  there 
exist  an  integer  N  =  N{e)  and  constants  a  , aN  such 
that 

N 

11  <P  -  l  an<Pjl  <  e  . 
n=o 

In  an  inner  product  space  X,  which  can  be  regarded  as 
normed  space  with  norm 


a  closed  sequence  is  complete,  but  the  reverse  is  not  true 
in  general. 


In  a  Hilbert  space,  however,  the  properties  of  closure 
and  completeness  are  equivalent . 


In  the  special  case  that  {y>*}  is  an 

n  n=o , 1 , . . . 

orthonormal  system  in  a  Hilbert  space  X  the  following 
statements  are  equivalent: 
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(i)  {tp*>  is  an  orthonormal  Hilbert-basis 

n  n=o,  1 , . . . 

(ii)  {cp£}  is  close<^  in  x 

(iii)  {<p*}  is  comPlete  in  x 

(iv)  For  each  ip  6  X, 

oo 

Hcpl!2  =  l  !  (<P,<P*)  |2 
n=o 

is  valid. 

(v)  For  each  (p  €  X  the  Fourier  expansion 

OO 

V  ■  [  (<p,<p*)<p* 
n=o 

holds  (in  sense  of  the  topology  of  X) . 

(vi)  For  ip(!  6  X, 

OO 

(<p,'n  =  l  (<p,tp*)  (v,4j) 
n=o 

("  Parseval’s  identity  ") 

Proofs  of  these  equivalences  can  be  found  e.g.  in  Davis 
(1963)  . 


7.  Trial  Functions  /  Base  Systems 


Some  base  systems  particularly  important  in  geodetic 
applications  should  be  mentioned  now: 

(i)  solid  spherical  harmonics 


Let  (S_  •} 


n,J  j=1,...,2n+1 


denote  a  maximal  linearly 


independent  system  of  (Laplace's)  surface  spherical 

harmonics  of  Oder  m.  Let  {H  .}  be  the 

n=o ,  1  , .  . . 
j= 1 , . . . , 2n+1 

system  of  solid  spherical  harmonics  defined  by 


Hnij(x)  =  |x|-|n+1)snij(5) 


(X  =  j  X  [  =  1)  . 


Then  the  system 

{h  .Is} 


n-o /!/••• 
j-*\  , . . . ,  2n+1 


(7.1) 


is  linearly  independent  and  complete  in  the  Hilbert 
space  *2(S). 


Proof . 

The  linear  independence  is  a  consequence  of  the  linear 
independence  of  the  system  {Sn  ^ }  of  surface  spherical 
harmonics  .. 

Therefore  it  remains  to  prove  that  {H  .|S> 

n,:3  n=o,1  , . . . 

j=1 , . . . , 2n+ 

is  complete  (and  consequently  closed)  in  the  Hilbert 
space  ^  (S)  . 


Suppose  that  ip  is  a  square-integrable  function 
2 

(<p  €  (S))  satisfying 

(4>,Hn  .)  =  /  ip(x)  H  .  (x)dS  =  0 

H/D  3  nf] 

for  all  n  =  o,1,...  ;  j  =  1,...,  2n+1  . 

We  remember  that  the  series  expansion 


wu 

-1—  =  4tt  l  jl 
|  x_y j  n£0  2n 


n  2n+1 


2n+1  j  ,n+1 
n=o  | y I 


l  sn  . 

j=1  n ' 3 


(?)  snfj(n) 


(x  *  | x |  ?  ,  y  =  |y  |  n,  ?2  =  n2  =  1 ) 


is  analytic  in  the  domain  given  by 


0  <  lxl  <  r0  '  ro  <  po/po  =  inflyi 

y€S 


\  x  6  E±  ,  y  €  S  > 

For  all  x  €  Ei  with  ]x|  <  rQ  we  thus  find  by  virtue 
of  our  assumption 


U1  (x)  =  /  ip(y)  j~yj  dS(y) 

s 

00  -  2] 
*  47r  l  "5TTIT  Wn 


l  2HTT  >x!n  Z  S  .(?)/<p(y)H  .(y)dS(y) 

n=o  ^n+l  j=1  n'3  S  n'3 


Analytic  continuation  shows  that  the  single-layer 
potential  U., (x)  vanishes  at  each  point  x  €  E. . 
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In  other  words,  the  equations 


(x  -  tn (x) )  =  0 

3U 


3n  (x) 


(x  -  xn(x) )  =  0 


hold  for  all  x  €  S  and  all  sufficiently  small  t  >  0. 


But  this  yields  using  the  relations  (4.8)  and  (4.9) 


lim  /  |U- (x  +  xn(x) ) | 2dS (x)  =  0 

t-*-o  S 

t-*-o 

3Ui  9 

lim  /  I  a~7~T  (x  +  xn(x))  +  4ir<p  (x)  |  dS  (x)  =»  0 
T-*o  S  1  3n(x) 

T-^O 


and 


lim 
T  -*-0 
T-*-0 


/  I  U’  (x)  +  2tr<p  (x)  |  2dS  (x) 
S  1 


0  . 


The  last  equation  can  be  rewritten  in  the  explicit  form 


/ 

S 


/ 

s 


<p(y)  a'nTxT  |''x-y~r  dS(y)  +  2^(x)  !2dS(x) 


o . 


In  other  words  we  have 

■  sr  l  v(y>  izm  ds(y>  =  vtx> 

in  the  sense  of  *?f2{S). 
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Iteration  yields 
2 


'■  37  >  ll  ~V“  m 7T  1^7-T*!y,)*(y2' 

t?  O  *  \  I 


=  ip(x) 


and 


(~  2it)  SSS^(Y3)  3n(y2}  I 3n(yi}  ly^y^T  3n(x)  'iy^xf ) 


i  a  i 


=  <P(x) 


dS(y1)dS(y2)dS(y3) 


Induction  states  that 


n 


s 


^{4>{yn)3n(yn_7)  |yn-y~_7T  *  *  *' sCTT  Ty^T^  (y1 } 


...  dS(yn) 

=  tf(x) 

2 

cor  all  integers  n  (in  the  sense  of  £  (S) )  .  Since 
iteration  reduces  the  order  of  singularity  (cf.  Kellogg 
(1929))  the  left  hand  side  is  a  continuous  function  of 
variable  x  at  S  (n  ?  0$  Thus  the  function  cp  can  be 

replaced  by  a  function  ip  €  (S)  satisfying 

/|ip(x)  -  tp(x)  !  2  dS  =  0  , 

~S  _y  2 

i.e.  tp  *  ip  in  the  sense  of  X  (S)  . 


For  the  continuous  function  ip,  however,  the  classical 

m 

limit  and  jump  formulas  are  valid: 


lim  U ^  (x  +  in (x) )  =  0 

x-*-o 

T-*-0 


lim 

T-*-0 

T-*-0 


au., 

fiilx) 


(X 


tn  (x)  ) 


4 rrtp  (x)  . 
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The  uniqueness  theorem  of  the  exterior  Dirichlet 
problem  then  yields  U-j  (x)  =  0  identically  for  all 
x  €  Ee.  But  this  means  that  =  0.  Consequently,  43  =  0 
on  the  surface  S.  q 

(ii)  singularity  functions  (mass  points) 

Let  F  be  a  (regular)  surface  in  the  inner  region 
E^  with  dist(F,S)  >  A  >  0. 


Suppose  that  {xn}  is  a  dense  point- 

n — o ,  1  ,  * .  . 

sequence  in  F(i.e.:  to  every  e  >  0  and  every  point 
x  €  F  there  exists  a  point  xn  such  that  jx-x  |  <  e) . 


Denote  by 


>n(x)  *  T^^TT  <xn  £  F  =  Ei> 

'  n ' 


(712) 


the  singularity  function  (fundamental  solution) 
at  xn  €  F. 

Then  the  system 


UJS} 


n=o , 1  , .  .  . 


is  linearly  independent  and  complete  in  the 
Hilbert  space  (S)  . 


Proof . 

Provided  that  xR  *  if  n  *  m  we  are  immediately  able 

to  prove  the  linear  independence  of  the  system 

Un|S}  in  st2(S). 

n-o / 1 / • • • 

Our  purpose  is  to  verify  the  completeness  of 
Un|S}  in  A ?  2(S)  . 


To  this  end  we  again  consider  a  function  (S)  . 

We  require  that 

!  VP (x)  0  (x)dS  =  0 

s 

for  n  =  0, 1 , • . .  . 

The  single-layer  potential 

U1  (x)  =  /  <x>  (y)  i" jy  |  ds  (y) 

s 

then  vanishes  at  all  points  x  €  F. 

n 

As  U1  is  continuous  on  the  surface  F ,  the  density  of 

the  system  {x  }  in  F  implies 

n=o , 1 , . . . 

U1  (x)  =  0 

for  all  x  €  F.  The  same  arguments  as  given  before  then 
assure  that  U1  vanishes  identically  in  R3 .  But  this 
yields  ip  =  0  in  ^  (S).  O 
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(iii)  "multipoles" 


Let  xQ  be  a  fixed  point  in  R  -  E  .  Denote  by 

{ <j)  }  the  sequence  of  "multipoles" 

n  n=o , 1 , . . . 
defined  by 


(»f  1 

3xo  lx-xO'(al=n 


(7.3) 


where 


[  a]  =  a-|  +a2+C£3 


ol  a-, 

■  (^r>  =  <db>  <s^J  • 

°  12  3  It 


Then  the  system  (<j>  is}  is  linearly 

n=o f 1 f  » . .  2 

independent  and  complete  in  the  Hilbert  space  x'iS) 


Proof . 

It  is  clear  that  the  potential  U1  with  all  its 
derivatives  of  arbitrary  order  vanishes  in  xq.  For  an 
analytic  function  this  implies  the  vanishing  in  a 
neighbourhood  of  the  point  xQ,  hence  in  -  Eg.  The 
next  steps  are  exactly  those  presented  above.  CD 

Other  types  of  trial  functions  have  been  listed  in 
Freeden  -  Kersten  (1982). 


Remark;  From  theory  of  spherical  harmonics  (cf.  e.g. 
Freeden (1 979) )  it  is  well-known  that  every (surf ace) 
spherical  harmonic  Sn  may  be  expressed  in  the  form 


V«>  ’  aj  V'V 


(I5I-D  , 


where  Pn  is  the  Legendre-polynomial  of  degree  n  and 
n1 '  * '  •  •  n2n+1  is  a  systein  points  n.  with  |n.  [  =  1  for 


which 
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a®t(VW).  .  ,  „  *  °- 

i=1 , . . . , 2n+1 
£=1 , .  .  . ,2n+1 


According  to  Maxwell's  representation  theorem 


(nV) 


n 


lx 


n  P_(5n) 
-  {-1}  nl  TTn+T 


the  Legendre  polynomials  may  be  obtained  by  repeated 
differentiations  of  the  fundamental  solution  of  the 
Laplace  equation  in  the  direction  of  the  unit  vector  n  . 
But  this  yields  (cf.  Cl.  Muller  (1966)) 


sn(s) 
i ..  i  n+1 


(-1) 

n! 


n 


2n+1 

I 

j-1 


(n.7) 


n 


_1_ 

x 


In  other  words,  spherical  harmonics  may  be  described 
always  by  combinations  of  multipoles.  Therefore  any 
approximation  by  spherical  harmonics  may  be  interpreted 
as  approximation  by  combinations  of  multipoles. 


8 .  Non  -  constructive  Approximation  Theorem 

In  the  following  we  summarize  our  results  developed 
in  Chapter  7  and  give  some  first  applications; 

Let  us  denote  by 


one  of  the  following  linearly  independent  systems  of 
functions  <j>n  €  P(Ee)  : 

(i)  solid  spherical  harmonics 

|xf  (n+1)  s  (c)  ,  (x-|x|C,  C2=1)  , 

11  •  j 

where  S  .  represents  a  maximal  linearly 
n,  j 

independent  system  of  (Laplace's)  surface 
spherical  harmonics. 


(ii)  "mass  points"  (reciprocal  distances) 

_ 

|  x-x  r  '  Xn  ^  E  '  E  *  Ei  > 
i  n  1  n 


where  {x  ) 
n 

i  .  e.  :  (a) 
(8) 


is  fundamental  in  E, 

n  o /I/... 

dist(xn,S)  >  0  for  all  n 
for  each  <p  €  cf0)  (EuS)  n  C  ( 2)  (E  ) 


with  A0  *  0  in  E  the  assumption 
<Mxn)  =  0,n=o,1,...,  implies  <t>  I E  = 


(iii)  "multipoles" 
z  0L 
(^'o) 


X  -X, 


[  a]  =n 


9 


where  x^  is  a  point  of  the  complement  of 
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Then  the  sequence 


(<Pn> 

n  n=o , 1 


§  •  •  • 


of  restrictions 

<Pn  =  ^nl  s  '  n=° / 1 , * . .  , 

is  a  closed  system  in  the  Hilbert  space  2 ( S)  ,  i.e.  to 
every  <p  €  if2(S)  and  every  e  >  o  there  exist  an  integer 
N=N(e)  and  coefficients  a  . . . . , a„  such  that 

O  N 

/  n  " 

J  /  l«P(x)  -  l  an<pn(x)  |MS  <  e. 

V  S  n-o  1 


In  particular,  to  every  <p  €  C(o)  (S)  and  every  e  >  0,  there 
exist  an  integer  N  =  N(e)  and  coefficients  aQ,...,aN  such 
that 


/ 


i  <P(x) 


N 

“  I 

n=o 


V>n(x) 


'dS  <  e 


(Observe  the  space  C^o) (S)  is  subset  of  the 
space  Jc  (S)  )  . 


By  application  of  the  regularization  theorem  we  obtain: 

Let  <j>  be  the  (unique)  function  of  class  P(E  )  with 
4>  |  S  *  '  (S)  •  Then,  to  every  e  >  0,  there  exist 

an  integer  N  =  N(e)  and  coefficients  a  ,  ...,aN  such  that 
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N 

sug  |  <J>  (x)  -  l  an<j>n(x)  | 
x€B  n=o 


/ 


j  k(y) 

s 


N 

-  I 

n=o 


am„  (y) 
n  n 


'dS(y) 


<  C  *  e 

for  each  compact  subset  B  c  Eg  with  dist  (B,S)  >  X  >  0. 

Unfortunately,  the  theorem  is  non  -  constructive, 
because  we  have  no  further  information  about  the 
coefficients  aQ,...,aN-  In  order  to  derive  constructive 
approximation  theorems  we  have  to  orthonormalize  the 
system  {(p  }  by  virtue  of  the  well-known 

Gram-Schmid?=<?6rtfi6norrnalization)  process . 
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9.  Constructive  Approximation  Theorem 


Let  {d»  }  „  be  one  of  the  systems  of  harmonic 

n  n=o, 1 , . . . 

functions  <J>  €  P(E  )  listed  in  Chapter  8  (i.e.  solid 

n  e 

spherical  harmonics/  mass  points,  multipoles) .  Let  us 
denote  by  as  usual  the  sequence  of 


n=o , 1 , . . . 


restrictions  tpn  =  <^nl  S,n=o,1 , . . .  .  Then,  corresponding  to 

the  countably  infinite  sequence  (ip  }  a  system 

n  n=o , 1 , . . . 

{<p*}  of  functions  <p*  satisfying 

n  n*o , 1 , . . .  n 


!  w*  (y)<p£(y)dS  (y) 

s 


(9.1) 


can  be  determined  recursively  in  the  following  way 
(cf.  Davis  (1961)  Sec. 8. 3) 


irr  '  *0  *  *0  '  ^o*2  =  (W 


*1  m  WJ  '  *1  -  *1  -  ^I'^o 


^  - 


ITT  '  ^2  ~  ^2  (^2,tpk)lpk 

£  K— O 


03* 

n+1 


*n+1  '  »n+1  ‘VrtK 


("Gram  -  Schmidt 


orthonoirmalization" ) 


According  to  the  exterior  Dirichlet  problem  of  potential 
theory,  there  exists  one  and  only  one  function  <j>^  6  &  (Ee) 

satisfying  <j>*|S  =  <$*,  n  =  o,1,...  . 

Therefore,  as  application  of  Corollary  2  (Chapt.  5^  we 
obtain: 


Let  4>  be  the  (unique)  function  of  class  ^(E  ) 

with  <p  (  S  =  cp,cp  €  C(o)  (S)  .  Then  the  £  2 -convergence  of 
the  sum 


implies  ordinary  pointwise  convergence  of  the  sum 


N  Nr 

l  (<p,<p*) <J>*  (x)  =  l  I  4>(y)  (y)ds  4>*(x) 
n-o  n=oLS  -I 


to.  4>(x)  as  N  -*■  00  for  each  x  E  Eg.  On  each  compact 
subset  B  with  B  c  Ee  and  dist(B,S)  >  X  >  0  the 
convergence  is  uniform: 

N 

l  (4>,M>*)  4>*  (x)  -  4>(x) 
n=o 


lim  sup 
n-»°°  xeff 


0 


The  approximation  technique  proposed  here,  indeed,  is 
constructed  so  as  to  have  the  so  -  called  permanence 
property: 

The  transition  from  the  N-th  truncation 
N 

I  <t>*  (x) 

n=o 

to  (N+1 ) -th  sum 

N+1 

l  (x) 

n=o 

necessitates  merely  the  addition  of  one  more  term, 
all  the  terms  obtained  formerly  remaining  unchanged. 

The  price  to  be  paid  for  the  convenience  of  the 
permanence  property  is  the  process  of  orthonormalization. 


The  theorem  is  constructive,  since  the  numbers 


(cpf(p*)  =  /  <p(y)<p*  (y)dS  (y) 

s 

are  the  Fourier  coefficients  of  with  respect  to  the 
system  {<p*} 

n=o, 1 , . . .  * 


(9.2) 


It  is  clear  from  the  structure  of  the  recursion  process 
that  ^n+1 »  and  hence,  <p*+1  ,  is  a  linear  combination  of 

<Pn+1(cf.(9  .1)).  in  order  to  determine  the  linear 
combination  tp*+1  from  tpQ,  —  ,<f>n+1  we  have  to  evaluate 
scalar  products  and  norms,  i.e.  surface  integrals  over  S. 
Assumed  the  surface  S  is  available  in  explicit  form  (e.g. 
in  the  case  of  sphere,  ellipsoid,  spheroid)  a  direct  attack 
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to  the  integrals  occuring  in  the  orthonormalization 
procedure  seems  possible  using  efficient  methods  of 
numerical  quadrature  or  cubature. 

Furthermore,  we  have 


—  N 

/U<y)|  ds  -  l  /  4>  (y)  <f>*  (y)ds  [ 

c  n=o  S  n  1 


+ 


N 

l 


n=o 


/  <P  (y)  4>*  (y)  ds  ]  2 
s  n 


for  any  selection  of  constants  aQ,...,an. 


Since  the  first  two  terms  on  the  right  hand  side  are 
independent  of  the  coefficients  a  , ...,an,  the  problem 
of  finding  the  best  approximation  by  linear  combination 
of  <J)q  : 

N  - 

min  /|<J>(y)  -  [  a  4>(y)|  dS  (9.4) 

an  S  n=o 

is  solved  if  and  only  if 


an  =  /  <P  (y)<p*(y)ds. 


(9.5) 


In  other  words. 


N  „ 

/ U  (y)  -  l  (d>,<l)*)<l>*(y)  I  ds 

S  n=o 

N  , 

<  fl<P(y)  -  l  ai.*(y)  TdS 
S  n=o  “ 


(9.6) 
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for  any  selection  of  constants  a  . 

o  n 

This  is  characteristic  for  (Fourier)  expansions. 


49 


t- 


(iv)  unions  of  the  base  systems  mentioned  under 
(i)  -  (iii) 

2.  Choose  an  appropriate  "approximation"  to  the  earth's 
surface  S  (e.g.  sphere,  ellipsoid , spheroid ,  ...) 

3.  Orthonormalize  the  system 

{ <j>n)n=0, 1 , .  .  .  ,  ^n'^=ipn  ’ 

with  respect  to  the  boundary  S.  This  yields  a  closed 
and  complete  system 


* 

{<p  }n=0, 1  .  .  . 
n 

in  the  Hilbert  space  ^  (S) . 
Define  the  system 


* 

{ (J>  }  / 1  f  •  •  • 

n 

by  **  f.P( Ee)  ,  4>n !  S=tPn  . 

4.  Calculate  the  coefficients 

(V,<p*)  =  /  V(y)  (y)dS(v) 

il  ^  I  • 

by  numerical  integration  methods. 


r  .  ' 

b. 

f* 


5o 


5 .  the  sum 

N  r 

V(N)(x)  =  l  [f  V(y)cp*(y)dF(v)J  **<x) 

serves  as  (globally  valid )  anproximation  for  the 
gravitational  potential 

V  £  P(F  ) 

e 


Therefore  the  whole  approximation  process  essentially 
amounts  to  computing  various  types  of  surface  integrals 
ove  r  S . 


Remark:  For  a  spherical  earth’s  model 
E  =  { x€R3  |  !  x  I  <  ?.} 


(R:  mean  radius  of  the  earth)  and  the  base  system  of 
spherical  harmonics  our  (numerical)  procedure  coincides 
with  the  conventional  computation  of  the  !?-th  truncated 
series  expansion  of  a  harmonic  function  into  solid  spheri¬ 
cal  harmonics  (cf.  e.g.  Heiskanen-Moritz ,  Chapter  1.16, 
(1967) ) . 


From  a  geodetic  point  of  view  it  should  he  noted  that 
the  external  gravitational  potential  ^  of  the  earth  may  be 
uniformly  approximated  by  the  sum 


n  *  *  :i  r 

(X)  -  l  (7,<fl  )  <4  (x)  =  l  /  ' 

p_— O  11  n  n=D  L  c 


:t 

i 

n=0 


,7  (y)'Pn  (y)  d.c  (y)  J 
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in  the  sense  that  for  any  civen  e  >  o  thorp  exists  an  integer 
N  =  N (c )  such  that  the  relation 

^  P  r  *  n  * 

!  I  I  /  V(y)ij>  (y)dfi(y)<j>  (x)  -  V(x)  e 
n=0  S  J 

holds  uniformly  everywhere  outside  and  on  any  closed  surface 
completely  surrounding  the  (earth's)  surface  F  in  the  outer 
space  E  .  The  surrounding  surface  may  be  arbitrarily  close  tc 
the  earth's  surface. 

Consequently  the  external  gravitational  potential  XT  of  the 
earth  can  be  expanded  into  the  series 

l  fi  V  (y)ip*  (y)d£  (y)l  <j>  (y) 
n=0  L  S  n  J 

actually  outside  and  on  any  closed  surface  completely 
surrounding  the  earth's  surface  S  in  the  outer  space  E  .  In 
other  words,  we  have  guaranteed  the  possibility  of  an  ap¬ 
proximation  of  the  external  gravitational  potential  in  the 

whole  outer  soace  E  bv  a  constructive  version  of  ,uro<? ' s 

- -  -  e  -  ■  —  _  —  - - - 


theorem. 


Numerical  Approximation 


From  the  structure  of  the  Gram-Schmidt  orthonorma Sizing 

* 

process  we  know  that  the  functions  <>  are  related  to  the 

n 

functions  <()  by  the  following  linear  equations 


>  =  I  D  .  <j>  .  , 

n  jio  n,j  ’3 


(11.1 


(n=0, . . . ,M) 


•where  the  matrix 


D  ^Dn,j^  n=o,...,N 
j  =0  (•.•(*! 


(11.2 


is  a  lower  triangular  matrix  (cf .  Chapter  9)  .  Thus  it  is 
clear  that  there  are  several  ways  to  express  the  best 
approximation  $  ^  to  an  element  $  €  ^(Eg)  by  a  combi¬ 
nation  of  given  (independent)  trial  functions: 


(i)  as  a  linear  combination 


<X>  =J0an  Vx> 


(11.3 


of  the  given  elements 


o 

t  %  f 


(ii)  as  a  linear  combination 


of  the  orthonormalized  elements 


a  ^ 


In  order  to  exploit  the  permanence  property  ortho- 
normalization  is  indispensable.  Orthonormalizaticn ,  however, 
is  equivalent  to  the  determination  of  the  matrix  0,  and 
computation  of  the  matrix  D  is  nothing  else  than  the 
inversion  of  a  well-def ined  (triangular)  matrix  C 
(of.  (1 1 . 1 8) ) .  From  a  numerical  point  of  view,  however, 
inversion  is  often  not  very  economical. 

The  situation  is  different  if  we  are  interested  only  in 

(\j)  • 

one  best  approximation  <j>  of  special  (prescribed)  degree 
N.  In  this  case,  a  numerical  technique  can  be  used  without 
any  need  of  inverting.  This  technique  is  based  on  Cholesky's 
factorization  theorem  as  will  be  described  now. 

To  this  end  we  begin  with  the  proof  of  the  following 
theorem: 


Let  the  sum  given  by 
i  (N)  (x)  =  l  an  <t>n  (x) 


(11.5) 


be  the  best  approximation  to  *  from  amonc  the 
linear  combination  of  Then,  the 

O  ,  .  .  •  ,  i! 

coefficients  a_  a  are  the  solutions  of  tKe 

■  /  f  m  m  •  f  n 

following  linear  system 


(*0 '*0)  a0  + 


+  n  ^  av  ~ 


(0 


(1 1 .«) 


("}.,» <tv)  av  =  (6 , 6,7)  . 


ao 


("normal  equations") 
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Proof:  The  coefficients  a„  are  determined 

-  0, . . . ,  N 


such  that 


N  N  *  * 

l  <j>n(x)  =  l  (<(>  ,<i>n)  <}>n. 


_  n  n 
n=0  n=0 


(11.7) 


Now,  it  is  easy  to  see  that 


N  +  +  * 

(*”  £  ( 4>  ^  4>  ^  ^  r  4>  ) 

j  =0  3  3  n 


(11.8) 


-  ( 4>  /  4>  )  ~  ( 1 4>  ) 

n  n 


■  0 


for  n=0, . . . ,N . 


But  this  means  that 

N 

(<t>“l3i.(t|i,./<t>)=0  (11.9) 

k=0  K  K  n 

for  n=0, . . . ,N. 

In  other  words,  $  minus  its  best  approximation  by  linear 
combinations  of  4>_  <f>.7  is  orthogonal  to  each  a.  □ 

The  normal  equations  (11.6)  can  be  rewritten  in  the 
vectorial  form 

T 

Ga  =  b  ,  b  =  (bQ, . . • , bN)  (11.10) 


where  b  €  R  ,  is  the  (N+1)-tupel 


b 


I !4'V  \ 

\  J 


and  G  is  the  Gram  matrix 


j  ( W 


(*p' V\ 

(W 


(11.11) 


(11.12) 


As  symmetric  and  positive  definite  matrix,  G  can  be 
decomposed  uniquely  in  the  form 

G  =  C  -  CT  ,  (11.13) 

where  C  is  a  lower  triangular  matrix  with  positive  diagonal 
elements.  This  splitting  of  G  is  known  as  the  Choleskv 
factorization. 

~ “ M+1  T 

Hence,  the  solution  a  £  R  ,  a  =  (an  aM) ,  can  be 

found  by  the  following  direct  method: 

(i)  factorization  of  the  Gram  matrix: 

T 

G  =  C-C 


(ii)  computation  of  the  vector 


of  the  linear  system 

C  d  =  b 


by  forward  substitution 


(iii)  computation  of  the  vector 


a  €  RN+1  ,  aT  =  (a0 . a,,)  , 


of  the  linear  system 


CTa  =  d 


by  backward  substitution. 


Remark:  As  is  well-known  from  numerical  analysis,  the 

procedure  is  very  economical  (  as  regards  to  the 
algorithmic  operations)  and  very  stable.  The 
solution  a  6  RN+”'  is  obtainable  without  any  need 
of  inverting. 


Observing  the  representation  (11.1)  we  get 
N  N 


U 


*  * 


k'V  -  J0  Dk,lD4.j  (W  ■  Sk,/ 


(11.14) 


This  is  eauivalent  to 


DGDT  =  I 


(11.15) 


(I:  unit  matrix) 


Consequently , 


G  =  d“1-(D_1)T  .  (11.16) 

On  the  other  hand,  we  have  the  (unique)  decomposition 

G  =  C  •  CT  .  (11.17) 

Thus  a  comparison  between  the  Gram-Schmidt  orthonorr.alizing 
and  the  Cholesky  factorizing  process  shows  that 

C  -  D"1  .  (11.18) 

N  +  1  t 

Moreover,  the  vector  c  €  R  ,  c  =  (cn,...,cM)  ,  is 

-N  + 1  T  {  ^ 

equal  to  the  vector  d  €  R  ,  d'  =  (dn,...,dN): 


C,  *5  (  0  ,  <t  . )  —  d  .  • 

3  D  3 


(11.1Q) 


This  is  easily  seen  from  (11.1)  by  the  following 
implications: 


.-I 


Db  =  c  **  b  =  D  c 

*»  b  *  Cc 

*»  b  »  Cd 


(of.  ( ii) ) 


i.  e.  : 


c  =  d 


(11.20) 


Hence,  by  Cholesky' s  method,  the  best  approximation 

$  to  an  element  <J>  with  respect  to  the  (orthonormalized) 
*  *  _ 
system  <j>n  <j>.T  can  be  determined  without  explicit 

computation  of  the  system  { * ^ >  itself. 


In  order  to  get  an  impression  of  the  least  square 
method  proposed  here  we  discuss  some  simple  numerical 
examples : 


Our  aim  is  to  approximate  <p  by  best  approximations 

N 

(x)  =  l  a  4>  (x)  (12.1) 

k=1  K  K 

with  respect  to  the  system  of  (buried)  mass  points 


*k(x)  =  lx-x 


xk  €  Ei  . 


(12.2) 


(i)  Boundaries 

The  surface  S  we  choose  for  our  computations  are 
given  by  _ 

V  ^  v  ^ 

*1  A  ^  A  ^ 

S  =  {x  1  (— -)  +  (~)  +  (-f-)  =  1  1 

a  DC  ' 


with 


(a)  a  =  b  =  c  =  1 


(Example  1) 


(8)  a  =  1,25  ;  b  =  1  ;  c  =  0,75 

(Examples  2,3) 


By  the  equations 


sin  #  cos  ip  (0<tp  <  2tt  ) 

sin  &  sin  ip  (0  <  <  u  ) 


x.  =  a 


x2  =  b  sin  &  sin  ip 
=  c  cos  t9* 


the  rectangular  coordinates  x.j,X2,X2  are  related  to 
the  ellipsoidal  coordinates:  a, b,c  (  axes  ),  (polar 
distance)  ,  cp  (geocentric  longitude)  . 
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(iii)  Trial  functions 


As  trial  functions  we  use  the  singularity  functions 

1 


=  IX-*., 


*k  €  Ei. 


The  set  of  mass  points 


(12.3) 


is  generated  by  radial  projection  of  the  points  £, 
o  K 

N 


of  the  system  J*, 


-  {?1  '  *  •  ’  ,5N} 


(12.4) 


to  a  surface  S  .  The  systems  ^  are  unions  of  regular 
polyhedra: 

J* ^  -  {PS1}  (Tetrahedron) 

/g  =  {PS1}  V  {PS2 }  (Cube) 

^2o  =  ^®1)  w  {PS2}  \j  {PS3}  (Dodecahedron) 

■^32  =  ^2o  w  {PS4}  (Icosahedron/Dodecahedron) 

■^55  =  ^2  2  w  {PS5}  u  {PS6}  ( Icosahedra/Dodecahedra) 

where  the  point  systems  (PSi) ,  i=1,...,6,  are  given 
as  follows: 


Integration  method 


The  integrals  over  S  occuring  in  the  Gram  matrix  G 
were  computed  in  the  following  way: 

The  integrals  over  S  are  transformed  using  polar 
coordinates.  Then  all  the  computations  of  the 
(transformed)  integrals  are  performed  by  (iterated) 
application  of  Romberg's  quadrature  formula  (Romberg 
extrapolation).  Using  the  subroutine  DCADRE  (cf.  de 
Boor  (1977))  we  are  able  to  compute  both  the  approxi¬ 
mation  of  the  integral  and  the  estimation  of  the 
error  between  integral  and  quadrature  formula.  In  all 
the  cases,  the  actual  resulting  error  was  about  lo-  . 


Remark : 


For  the  Cholesky  factorization,  the  forward  and 
backward  substitution  powerful  routines  are  available 
(cf.  Dongarra  et.al . (1 979) ) . 
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(v)  error  estimates 


The  accuracy  of  the  best  approximation  $  ^  to  <j> 
will  be  investigated  by  discussing  the  errors 


LS  : 


MAR  : 


f  [({>( y)  -  4>(N)  (y)  ]  2  dS 
S 


(least  square  error) 


max 

x€S 


(x)  -  (ji  ^  (x) 
<}>  (x) 


(maximal  absolute 
relative  error) 

r  i~iW.  (y)  |  ds 

g  1  4>  (y)  l  as  . 


(mean  absolute  relative  error) 


The  error  LS  is  determined  by  the  Romberg 
integration  mentioned  above. 


The  errors  MAR,M  are  evaluated  for  a  discrete 
lattice  of  181  *  91  points  in  the  (q>  ,9>)  -plane  ( ip,  & 
polar  coordinates) . 


M  = 


1 

THTTTT 


181  91 

I  I 

i=l  J=1 


M 


#*ij> 


1} 


Example  1 : 


S  =  (x|x2  +  x2  +  x2 


~  /  1  2  2  2 

S  =  ( x | x  ^  +  x2  +  x3 


Boundary  values: 


(unit  sphere) 
7/10) 


<p(x)  =  I  x-y  |  ~  ,  y  =  -r  (1,1,1)  ,  x  e  S 


N 

LS 

MAR 

M 

4 

8 

0,25 

0,12  *  10-2 

0,80 

0,93  *  10’1 

0,24 

0,16  *  10”1 

2o 

0,50  •  10”3 

0,62  *  10"1 

0,77  ‘  10-2 

32 

0,23  •  10"3 

0,48  ’  10_1 

0,50  *  10-2 

56 

.3 

0,21  *  10 

0,46  *  10’1 

0,40  •  10-2 

2 


2 


1} 


Example  2: 


s  =  ( x  I  ( i r  25)  +  x2  +  (0,75} 


(ellipsoid) 

's'  -  { x  j x  =  rw  ,  w  €  S  ,  r  =  0 , 7  5 } 
Boundary  values: 

4>(x)  =  ’  x-y  1  ” 1  ,  y  -  (1,1,1)  ,  x  e  E 


N 

LS 

MAR 

M 

■ 

0,32 

1,12 

0,25 

8 

0,33* 10-2 

0,13 

o 

** 

to 

• 

—a 

o 

i 

-5 

-1 

2o 

O 

• 

o 

0,10 

0,12  •  10 

32 

O 

_ X 

• 

o 

1 

IsJ 

0,82*10-1 

o 

-J 

00 

O 

1 

to 

56 

0,99'10~3 

o 

00 
— * 

* 

o 

1 

0,59  •  10~2 

Example  3: 


S  =  {x'  (1  ,255  +  X2  +  Boris'1  ~  1  } 

(ellipsoid) 

S  *  (x  =  rw  ,  w  6  S  ,  r  =  0,75} 


Boundary  values: 


N 

LS 

MAR 

M 

■ 

1  ,37 

1  ,01 

0,42 

8 

1/25 

0,85 

0,42 

2o 

0,29 

0,42 

0;  1  3 

32 

0,53 • 10"1 

0,24 

o 

CJ 

00 

• 

o 

— 1 

5o 

O 

o 

• 

o 

1 

0,10 

0,19  *  10’1 
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0,98- 10-3 

0,35*10_1 

0,88  *  10"2 

Figure  6a:  error  cj>(x)  -  $  (x),  xeS  (N=4,  Example  2) 


$(x) 


xeS  (N**8,  Example  2) 


Figure  12a:  error  <(>(x)  -  $  (x) ,  xeS  (N=8,  Example  3) 


Figure  14a:  error  $(x) 


(x),  xeS  (N=32,  Example  3) 


xeS  (N=32 


Figure  15a:  error  <(>(x)  -  $  (x),  xeS  (N=56,  Example  3) 


(N) 

Figure  15b:  relative  error  <Kx)  ~  <|»  (x)  , 


xeS  (N«56,  Example  3) 


Final  Remarks : 


Mass  model  representations  have  been  shown  to  be 
adequate  tools  for  the  purpose  of  representing  the  external 
gravitational  potential.  The  relations  of  the  models  to 
the  physical  reajitiy  are  transparent:  the  harmonicity  of 
the  approximation  by  mass  model  representations  is 
guaranteed.  Point  masses  (reciprocal  distances)  are  easy 
to  handle,  reciprocal  distances  and  their  derivatives 
provide  simple  expressions  of  any  desired  gravitational 
field  quantity.  The  relation  between  point  masses  and 
spherical  harmonic  coefficients  of  the  external  potential 
is  straightforward  (cf .  Chapt.  7) .  By  combination  of  the 
base  functions  (e.g.  low  degree  spherical  harmonics  and 
mass  points)  physical  meanings  (center  of  mass,  moments 
of  inertia  etc.)  can  be  easily  implemented  into  the  model. 
Least  square  approximations  by  mass  model  representations 
can  be  deduced  for  not-necessarily  spherical  earth's 
models,  the  approximations  are  best  in  the  sense  of  the 
root  -  mean-square  error. 

There  are  of  course  some  drawbacks  in  the  approximation 
of  this  technique : theoretically  there  exist  infinite 
numbers  of  mass  distributions  compatible  physically  with 
the  earth's  gravitational  field?  in  practical 
applications,  however,  we  have  to  select  a  finite  number 
of  (multi) poles  which  are  both  economical  computationally 
and  relevant  geophysically.  Approximation  of  boundary 
values  and  external  gravitational  potential  <j»  is  achieved 
by  sums  <i>  ^  of  oscillating  character.  As  illustrated  in  our 
Example  3  (Chapter  12),  the  oscillations  often  grow  in 
number,  but  they  decrease  in  size  with  increasing  N. 
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Least  square  approximation  as  described  in  this  report 
provides  a  suitable  procedure  for  macro  modeling  by 
successive  oscillations.  It  is  not  a  technique  of 
osculating  character  (as,  for  example,  the  Taylor 
series) . 


Remark  :  A  detailed  discussion  of  oscillating  and  osculating 
approximation  in  classical  Fourier  theory  has  been  given  by 
A.  Sommer f eld  (1947) :  PARTIELLE  DIFFERENT IALGLEICHUNGEN  DER 
PHYSIK,  Vorlesungen  tiber  Theoretische  Physik,  Bd.  6,  Akademische 
Verlagsgesellschaft,  Leipzig  (6.Auflage  (1966)). 
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